R/Clink et al Variance Components Model and Figure Code.R
In DenaJGibbon/geographic_tarsiers:
# Tarsier MANOVA
### Load required libraries
library(corrplot)
library(bayesplot)
library(rstan)
rstan_options(auto_write=TRUE)
options(mc.cores = parallel::detectCores())
library(stringr)
library(viridis)
library(ggplot2)
library(stringr)
# Read in the data
tarsier.feature.table.f <- read.csv("/Users/denasmacbook/Geographic tarsiers/tarsier.female.march28.csv")
# Isolate features of interest
d.manova <- tarsier.feature.table.f[,c("note.1.dur", "note.2.dur","note.1.max.freq",
"note.2.max.freq","term.f","n.notes","note.rate")]
## Check the structure of the data
str(d.manova)
pairs(log(d.manova))
cor(d.manova)
## Log transform data
d.manova <- log(d.manova)
hist(d.manova$note.1.dur)
### Set-up data to pass to Stan.
# Integer-coded vector of group IDs
group.int <- as.numeric(tarsier.feature.table.f$group)
# Check structure
table(group.int)
# Integer-coded vector of site IDs
site.int <- as.numeric(as.factor(tarsier.feature.table.f$site))
# Check structure
table(site.int)
# Center data matrix at feature means
col.means <- apply(d.manova, MARGIN=2, FUN="mean")
y.centered <- sweep(d.manova, MARGIN=2, STATS=col.means)
# Create a data list to pass to Stan
data_list <- list(
K = dim(d.manova)[2],
J= length(unique(tarsier.feature.table.f$group)),
M= length(unique(tarsier.feature.table.f$site)),
N= dim(d.manova)[1],
y= as.matrix(y.centered), ## features centered at zero
group= group.int,
site= site.int
)
# Code to run the STAN mdoel
# NOTE: the .stan file must be linked in the code below
m.stan.tarsier.f.seven.features = stan(file="/Users/denasmacbook/Geographic tarsiers/tarsier.stan",
model_name = "m.stan.tarsier.f",
data=data_list, iter=1500, warmup=1000, chains=2,
cores=2,
control = list(stepsize = 0.5, adapt_delta = 0.99, max_treedepth = 20))
#save(m.stan.tarsier.f, file = "/Users/denasmacbook/Desktop/m.stan.tarsier.f.rda")
#load("/Users/denasmacbook/Desktop/m.stan.tarsier.f.rda")
## Check model output
# Create traceplots to check for mixing for site level variance
draws <- as.array(m.stan.tarsier.f.seven.features, pars="ICC_site")
mcmc_trace(draws)
stan_dens(m.stan.tarsier.f.seven.features, pars=c("ICC_site"))
round(summary(m.stan.tarsier.f.seven.features, pars=c("ICC_site"))$summary, 3)
# Create traceplots to check for mixing for group level variance
draws <- as.array(m.stan.tarsier.f.seven.features, pars="ICC_group")
mcmc_trace(draws)
stan_dens(m.stan.tarsier.f.seven.features, pars=c("ICC_group"))
round(summary(m.stan.tarsier.f.seven.features, pars=c("ICC_group"))$summary, 3)
# Check degrees of freedom parameter
stan_dens(m.stan.tarsier.f.seven.features, pars=c("DF_obs"))
round(summary(m.stan.tarsier.f.seven.features, pars=c("DF_obs"))$summary, 3)
stan_dens(m.stan.tarsier.f, pars=c("DF_group"))
round(summary(m.stan.tarsier.f, pars=c("DF_group"))$summary, 3)
stan_dens(m.stan.tarsier.f, pars=c("DF_site"))
round(summary(m.stan.tarsier.f, pars=c("DF_site"))$summary, 3)
post.samples <- rstan::extract(m.stan.tarsier.f.seven.features,
pars= c("ICC_site", "ICC_group", "Maha_sqd", "DF_site","DF_group",
"DF_obs", "Vcov_site", "Vcov_group", "Vcov_obs"),
permuted=TRUE)
# Create vector with feature names
features <- c("note.1.dur", "note.2.dur","note.1.max.freq",
"note.2.max.freq","term.f","n.notes","note.rate")
# Add column names to samples
colnames(post.samples$ICC_group) <- features
colnames(post.samples$ICC_site) <- features
# Convert to dataframe
icc.site <- as.data.frame(post.samples$ICC_site)
icc.group <- as.data.frame(post.samples$ICC_group)
# Calculate ICC for observation-level and add column names
icc.obs <- as.data.frame(1 - icc.group - icc.site)
colnames(icc.obs) <- features
### Code to make Figure 4. Posterior densities for the intraclass correlation coefficients for the three levels in our dataset
### (call, fefemale and site) for each feature of the tarsier duet phrase
## Note 2 duration
note.1.dur.site <- cbind.data.frame(icc.site$note.1.dur, rep("Site",length(icc.site$note.1.dur)))
note.1.dur.group <- cbind.data.frame(icc.group$note.1.dur,rep("Group",length(icc.group$note.1.dur)))
note.1.dur.obs <- cbind.data.frame(icc.obs$note.1.dur,rep("Obs",length(icc.obs$note.1.dur)))
colnames(note.1.dur.site) <- c("note.1.dur","icc")
colnames(note.1.dur.group) <- c("note.1.dur","icc")
colnames(note.1.dur.obs) <- c("note.1.dur","icc")
note.1.dur.dens.df <- rbind.data.frame(note.1.dur.site,note.1.dur.group,note.1.dur.obs)
head(note.1.dur.dens.df)
female.note.1.dur.plot <- ggplot(note.1.dur.dens.df, aes(x=note.1.dur, fill=icc)) + geom_density()+xlim(0,1)+#+ylim(0,4)+
scale_fill_manual(name = "icc",
values = c(alpha("blue", .5),
alpha("#6DCD59FF", .5),
alpha("#440154FF", .5)))+ theme_classic() +
guides(fill=F)+xlab("")+ylab("") +theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
axis.text.x = element_text(size=18))
female.note.1.dur.plot
## Note 2 duration
note.2.dur.site <- cbind.data.frame(icc.site$note.2.dur, rep("Site",length(icc.site$note.2.dur)))
note.2.dur.group <- cbind.data.frame(icc.group$note.2.dur,rep("Group",length(icc.group$note.2.dur)))
note.2.dur.obs <- cbind.data.frame(icc.obs$note.2.dur,rep("Obs",length(icc.obs$note.2.dur)))
colnames(note.2.dur.site) <- c("note.2.dur","icc")
colnames(note.2.dur.group) <- c("note.2.dur","icc")
colnames(note.2.dur.obs) <- c("note.2.dur","icc")
note.2.dur.dens.df <- rbind.data.frame(note.2.dur.site,note.2.dur.group,note.2.dur.obs)
head(note.2.dur.dens.df)
female.note.2.dur.plot <- ggplot(note.2.dur.dens.df, aes(x=note.2.dur, fill=icc)) + geom_density()+xlim(0,1)+#+ylim(0,4)+
scale_fill_manual(name = "icc",
values = c(alpha("blue", .5),
alpha("#6DCD59FF", .5),
alpha("#440154FF", .5)))+ theme_classic() +
guides(fill=F)+xlab("")+ylab("") +theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
axis.text.x = element_text(size=18))
female.note.2.dur.plot
## Note 1 Max Freq
note.1.max.freq.site <- cbind.data.frame(icc.site$note.1.max.freq, rep("Site",length(icc.site$note.1.max.freq)))
note.1.max.freq.group <- cbind.data.frame(icc.group$note.1.max.freq,rep("Group",length(icc.group$note.1.max.freq)))
note.1.max.freq.obs <- cbind.data.frame(icc.obs$note.1.max.freq,rep("Obs",length(icc.obs$note.1.max.freq)))
colnames(note.1.max.freq.site) <- c("note.1.max.freq.samples","icc")
colnames(note.1.max.freq.group) <- c("note.1.max.freq.samples","icc")
colnames(note.1.max.freq.obs) <- c("note.1.max.freq.samples","icc")
note.1.max.freq.dens.df <- rbind.data.frame(note.1.max.freq.site,note.1.max.freq.group,note.1.max.freq.obs)
head(note.1.max.freq.dens.df)
female.note.1.max.freq.plot <- ggplot(note.1.max.freq.dens.df, aes(x=note.1.max.freq.samples, fill=icc)) + geom_density()+xlim(0,1)+#+ylim(0,4)+
scale_fill_manual(name = "icc",
values = c(alpha("blue", .5),
alpha("#6DCD59FF", .5),
alpha("#440154FF", .5)))+ theme_classic() +
guides(fill=F)+xlab("")+ylab("") +theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
axis.text.x = element_text(size=18))
female.note.1.max.freq.plot
# Note 2 Max Freq
note.2.max.freq.site <- cbind.data.frame(icc.site$note.2.max.freq, rep("Site",length(icc.site$note.2.max.freq)))
note.2.max.freq.group <- cbind.data.frame(icc.group$note.2.max.freq,rep("Group",length(icc.group$note.2.max.freq)))
note.2.max.freq.obs <- cbind.data.frame(icc.obs$note.2.max.freq,rep("Obs",length(icc.obs$note.2.max.freq)))
colnames(note.2.max.freq.site) <- c("note.2.max.freq.samples","icc")
colnames(note.2.max.freq.group) <- c("note.2.max.freq.samples","icc")
colnames(note.2.max.freq.obs) <- c("note.2.max.freq.samples","icc")
note.2.max.freq.dens.df <- rbind.data.frame(note.2.max.freq.site,note.2.max.freq.group,note.2.max.freq.obs)
head(note.2.max.freq.dens.df)
female.note.2.max.freq.plot <- ggplot(note.2.max.freq.dens.df, aes(x=note.2.max.freq.samples, fill=icc)) + geom_density()+xlim(0,1)+ylim(0,4)+
scale_fill_manual(name = "icc",
values = c(alpha("blue", .5),
alpha("#6DCD59FF", .5),
alpha("#440154FF", .5)))+ theme_classic() +
guides(fill=F)+xlab("")+ylab("") +theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
axis.text.x = element_text(size=18))
female.note.2.max.freq.plot
## Note 1 max freq
term.f.site <- cbind.data.frame(icc.site$term.f, rep("Between-site",length(icc.site$term.f)))
term.f.group <- cbind.data.frame(icc.group$term.f,rep("Between-fefemale",length(icc.group$term.f)))
term.f.obs <- cbind.data.frame(icc.obs$term.f,rep("Within-fefemale",length(icc.group$term.f)))
colnames(term.f.site) <- c("term.f.samples","icc")
colnames(term.f.group) <- c("term.f.samples","icc")
colnames(term.f.obs) <- c("term.f.samples","icc")
term.f.dens.df <- rbind.data.frame(term.f.site,term.f.group,term.f.obs)
head(term.f.dens.df)
female.term.f.plot <- ggplot(term.f.dens.df, aes(x=term.f.samples, fill=icc)) + geom_density()+xlim(0,1)+ #ylim(0,90)+
scale_fill_manual(name = "ICC",
values = c(alpha("blue", .5),
alpha("#6DCD59FF", .5),
alpha("#440154FF", .5)))+ theme_classic() +
guides(fill=F)+xlab("")+ylab("") +theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
axis.text.x = element_text(size=18),legend.text=element_text(size=24))+
xlab("")+ylab("") + theme(legend.position = c(.8,.87))
female.term.f.plot
## Note 2 duration
n.notes.site <- cbind.data.frame(icc.site$n.notes, rep("Island",length(icc.site$n.notes)))
n.notes.group <- cbind.data.frame(icc.group$n.notes,rep("Female",length(icc.group$n.notes)))
n.notes.obs <- cbind.data.frame(icc.obs$n.notes,rep("Phrase",length(icc.group$n.notes)))
colnames(n.notes.site) <- c("n.notes.samples","icc")
colnames(n.notes.group) <- c("n.notes.samples","icc")
colnames(n.notes.obs) <- c("n.notes.samples","icc")
n.notes.dens.df <- rbind.data.frame(n.notes.site,n.notes.group,n.notes.obs)
head(n.notes.dens.df)
female.n.notes.plot <- ggplot(n.notes.dens.df, aes(x=n.notes.samples, fill=icc)) + geom_density()+xlim(0,1)+ #ylim(0,90)+
scale_fill_manual(name = "ICC",
values = c(alpha("blue", .5),
alpha("#6DCD59FF", .5),
alpha("#440154FF", .5)))+ theme_classic() +
guides(fill = guide_legend(size=10, keywidth = 3, title="",keyheight = 2))+ xlab("")+ylab("") +theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
axis.text.x = element_text(size=18))+ theme(legend.text=element_text(size=12,face = "bold"))
female.n.notes.plot
## Note 2 duration
note.rate.site <- cbind.data.frame(icc.site$note.rate, rep("Site",length(icc.site$note.rate)))
note.rate.group <- cbind.data.frame(icc.group$note.rate,rep("Group",length(icc.group$note.rate)))
note.rate.obs <- cbind.data.frame(icc.obs$note.rate,rep("Obs",length(icc.group$note.rate)))
colnames(note.rate.site) <- c("note.rate.samples","icc")
colnames(note.rate.group) <- c("note.rate.samples","icc")
colnames(note.rate.obs) <- c("note.rate.samples","icc")
note.rate.dens.df <- rbind.data.frame(note.rate.site,note.rate.group,note.rate.obs)
head(note.rate.dens.df)
female.note.rate.plot <- ggplot(note.rate.dens.df, aes(x=note.rate.samples, fill=icc)) + geom_density()+xlim(0,1)+ #ylim(0,90)+
scale_fill_manual(name = "ICC",
values = c(alpha("blue", .5),
alpha("#6DCD59FF", .5),
alpha("#440154FF", .5)))+ theme_classic() +
guides(fill=F)+xlab("")+ylab("") +theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
axis.text.x = element_text(size=18))
female.note.rate.plot
# Combine all plots
female.plot <- cowplot::plot_grid(female.note.1.dur.plot, female.note.2.dur.plot,female.note.1.max.freq.plot,female.note.2.max.freq.plot,
female.term.f.plot,female.n.notes.plot,female.note.rate.plot,
label_x=0.3,label_y = 1,
labels=c(" Note 1 Dur", "Note 2 Dur", "Note 1 Max Freq",
"Note 2 Max Freq","Terminal Freq","Number of Notes", " Note Rate"),
ncol=2)
######## Male MANOVA
male.tarsier.df <- read.csv("tarsier.feature.table.m.csv")
d.manova.male <- male.tarsier.df[,c("note.1.dur", "note.2.dur","note.1.max.freq",
"note.2.max.freq","term.f","n.notes","note.rate")]
## Check the structure of the data
str(d.manova.male)
pairs(d.manova.male)
cor(d.manova.male)
## Log transform data
d.manova.male <- log(d.manova.male)
### Set-up data to pass to Stan.
# Integer-coded vector of group IDs
group.int <- as.numeric(male.tarsier.df$group)
# Check structure
table(group.int)
# Integer-coded vector of site IDs
site.int <- as.numeric(as.factor(male.tarsier.df$site))
# Check structure
table(site.int)
# Center data matrix at feature means
col.means <- apply(d.manova.male, MARGIN=2, FUN="mean")
y.centered <- sweep(d.manova.male, MARGIN=2, STATS=col.means)
# Create a data list to pass to Stan
data_list <- list(
K = dim(d.manova.male)[2],
J= length(unique(male.tarsier.df$group)),
M= length(unique(male.tarsier.df$site)),
N= dim(d.manova.male)[1],
y= as.matrix(y.centered), ## features centered at zero
group= group.int,
site= site.int
)
# Code to run the STAN mdoel
# NOTE: the .stan file must be linked in the code below
m.stan.tarsier.m = stan(file="/Users/denasmacbook/Desktop/Clink et al Multivariate Analysis/MANOVA_MNG_oct18_2017.stan",
model_name = "m.stan.tarsier.m",
data=data_list, iter=1000, warmup=500, chains=2,
cores=2,
control = list(stepsize = 0.5, adapt_delta = 0.99, max_treedepth = 20))
## Check model output
# Create traceplots to check for mixing for site level variance
draws <- as.array(m.stan.tarsier.m, pars="ICC_site")
mcmc_trace(draws)
stan_dens(m.stan.tarsier.m, pars=c("ICC_site"))
round(summary(m.stan.tarsier.m, pars=c("ICC_site"))$summary, 3)
# Create traceplots to check for mixing for group level variance
draws <- as.array(m.stan.tarsier.m, pars="ICC_group")
mcmc_trace(draws)
stan_dens(m.stan.tarsier.m, pars=c("ICC_group"))
round(summary(m.stan.tarsier.m, pars=c("ICC_group"))$summary, 3)
# Check degrees of freedom parameter
stan_dens(m.stan.tarsier.m, pars=c("DF_obs"))
round(summary(m.stan.tarsier.m, pars=c("DF_obs"))$summary, 3)
stan_dens(m.stan.tarsier.m, pars=c("DF_group"))
round(summary(m.stan.tarsier.m, pars=c("DF_group"))$summary, 3)
stan_dens(m.stan.tarsier.m, pars=c("DF_site"))
round(summary(m.stan.tarsier.m, pars=c("DF_site"))$summary, 3)
post.samples <- extract(m.stan.tarsier.m,
pars= c("ICC_site", "ICC_group", "Maha_sqd", "DF_site","DF_group",
"DF_obs", "Vcov_site", "Vcov_group", "Vcov_obs"),
permuted=TRUE)
# Create vector with feature names
features <- c("note.1.dur", "note.2.dur","note.1.max.freq",
"note.2.max.freq","term.f","n.notes","note.rate")
# Add column names to samples
colnames(post.samples$ICC_group) <- features
colnames(post.samples$ICC_site) <- features
# Convert to dataframe
icc.site <- as.data.frame(post.samples$ICC_site)
icc.group <- as.data.frame(post.samples$ICC_group)
# Calculate ICC for observation-level and add column names
icc.obs <- as.data.frame(1 - icc.group - icc.site)
colnames(icc.obs) <- features
### Code to make Figure 4. Posterior densities for the intraclass correlation coefficients for the three levels in our dataset
### (call, female and site) for each feature of the tarsier duet phrase
## Note 2 duration
note.1.dur.site <- cbind.data.frame(icc.site$note.1.dur, rep("Site",length(icc.site$note.1.dur)))
note.1.dur.group <- cbind.data.frame(icc.group$note.1.dur,rep("Group",length(icc.group$note.1.dur)))
note.1.dur.obs <- cbind.data.frame(icc.obs$note.1.dur,rep("Obs",length(icc.obs$note.1.dur)))
colnames(note.1.dur.site) <- c("note.1.dur","icc")
colnames(note.1.dur.group) <- c("note.1.dur","icc")
colnames(note.1.dur.obs) <- c("note.1.dur","icc")
note.1.dur.dens.df <- rbind.data.frame(note.1.dur.site,note.1.dur.group,note.1.dur.obs)
head(note.1.dur.dens.df)
male.note.1.dur.plot <- ggplot(note.1.dur.dens.df, aes(x=note.1.dur, fill=icc)) + geom_density()+xlim(0,1)+#+ylim(0,4)+
scale_fill_manual(name = "icc",
values = c(alpha("blue", .5),
alpha("#6DCD59FF", .5),
alpha("#440154FF", .5)))+ theme_classic() +
guides(fill=F)+xlab("")+ylab("") +theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
axis.text.x = element_text(size=18))
male.note.1.dur.plot
## Note 2 duration
note.2.dur.site <- cbind.data.frame(icc.site$note.2.dur, rep("Site",length(icc.site$note.2.dur)))
note.2.dur.group <- cbind.data.frame(icc.group$note.2.dur,rep("Group",length(icc.group$note.2.dur)))
note.2.dur.obs <- cbind.data.frame(icc.obs$note.2.dur,rep("Obs",length(icc.obs$note.2.dur)))
colnames(note.2.dur.site) <- c("note.2.dur","icc")
colnames(note.2.dur.group) <- c("note.2.dur","icc")
colnames(note.2.dur.obs) <- c("note.2.dur","icc")
note.2.dur.dens.df <- rbind.data.frame(note.2.dur.site,note.2.dur.group,note.2.dur.obs)
head(note.2.dur.dens.df)
male.note.2.dur.plot <- ggplot(note.2.dur.dens.df, aes(x=note.2.dur, fill=icc)) + geom_density()+xlim(0,1)+#+ylim(0,4)+
scale_fill_manual(name = "icc",
values = c(alpha("blue", .5),
alpha("#6DCD59FF", .5),
alpha("#440154FF", .5)))+ theme_classic() +
guides(fill=F)+xlab("")+ylab("") +theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
axis.text.x = element_text(size=18))
male.note.2.dur.plot
## Note 1 Max Freq
note.1.max.freq.site <- cbind.data.frame(icc.site$note.1.max.freq, rep("Site",length(icc.site$note.1.max.freq)))
note.1.max.freq.group <- cbind.data.frame(icc.group$note.1.max.freq,rep("Group",length(icc.group$note.1.max.freq)))
note.1.max.freq.obs <- cbind.data.frame(icc.obs$note.1.max.freq,rep("Obs",length(icc.obs$note.1.max.freq)))
colnames(note.1.max.freq.site) <- c("note.1.max.freq.samples","icc")
colnames(note.1.max.freq.group) <- c("note.1.max.freq.samples","icc")
colnames(note.1.max.freq.obs) <- c("note.1.max.freq.samples","icc")
note.1.max.freq.dens.df <- rbind.data.frame(note.1.max.freq.site,note.1.max.freq.group,note.1.max.freq.obs)
head(note.1.max.freq.dens.df)
male.note.1.max.freq.plot <- ggplot(note.1.max.freq.dens.df, aes(x=note.1.max.freq.samples, fill=icc)) + geom_density()+xlim(0,1)+#+ylim(0,4)+
scale_fill_manual(name = "icc",
values = c(alpha("blue", .5),
alpha("#6DCD59FF", .5),
alpha("#440154FF", .5)))+ theme_classic() +
guides(fill=F)+xlab("")+ylab("") +theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
axis.text.x = element_text(size=18))
male.note.1.max.freq.plot
# Note 2 Max Freq
note.2.max.freq.site <- cbind.data.frame(icc.site$note.2.max.freq, rep("Site",length(icc.site$note.2.max.freq)))
note.2.max.freq.group <- cbind.data.frame(icc.group$note.2.max.freq,rep("Group",length(icc.group$note.2.max.freq)))
note.2.max.freq.obs <- cbind.data.frame(icc.obs$note.2.max.freq,rep("Obs",length(icc.obs$note.2.max.freq)))
colnames(note.2.max.freq.site) <- c("note.2.max.freq.samples","icc")
colnames(note.2.max.freq.group) <- c("note.2.max.freq.samples","icc")
colnames(note.2.max.freq.obs) <- c("note.2.max.freq.samples","icc")
note.2.max.freq.dens.df <- rbind.data.frame(note.2.max.freq.site,note.2.max.freq.group,note.2.max.freq.obs)
head(note.2.max.freq.dens.df)
male.note.2.max.freq.plot <- ggplot(note.2.max.freq.dens.df, aes(x=note.2.max.freq.samples, fill=icc)) + geom_density()+xlim(0,1)+#+ylim(0,4)+
scale_fill_manual(name = "icc",
values = c(alpha("blue", .5),
alpha("#6DCD59FF", .5),
alpha("#440154FF", .5)))+ theme_classic() +
guides(fill=F)+xlab("")+ylab("") +theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
axis.text.x = element_text(size=18))
male.note.2.max.freq.plot
## Note 1 max freq
term.f.site <- cbind.data.frame(icc.site$term.f, rep("Between-site",length(icc.site$term.f)))
term.f.group <- cbind.data.frame(icc.group$term.f,rep("Between-female",length(icc.group$term.f)))
term.f.obs <- cbind.data.frame(icc.obs$term.f,rep("Within-female",length(icc.group$term.f)))
colnames(term.f.site) <- c("term.f.samples","icc")
colnames(term.f.group) <- c("term.f.samples","icc")
colnames(term.f.obs) <- c("term.f.samples","icc")
term.f.dens.df <- rbind.data.frame(term.f.site,term.f.group,term.f.obs)
head(term.f.dens.df)
male.term.f.plot <- ggplot(term.f.dens.df, aes(x=term.f.samples, fill=icc)) + geom_density()+xlim(0,1)+ #ylim(0,90)+
scale_fill_manual(name = "ICC",
values = c(alpha("blue", .5),
alpha("#6DCD59FF", .5),
alpha("#440154FF", .5)))+ theme_classic() +
guides(fill=F)+xlab("")+ylab("") +theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
axis.text.x = element_text(size=18),legend.text=element_text(size=24))+
xlab("")+ylab("") + theme(legend.position = c(.8,.87))
male.term.f.plot
## Note 2 duration
n.notes.site <- cbind.data.frame(icc.site$n.notes, rep("Island",length(icc.site$n.notes)))
n.notes.group <- cbind.data.frame(icc.group$n.notes,rep("Male",length(icc.group$n.notes)))
n.notes.obs <- cbind.data.frame(icc.obs$n.notes,rep("Phrase",length(icc.group$n.notes)))
colnames(n.notes.site) <- c("n.notes.samples","icc")
colnames(n.notes.group) <- c("n.notes.samples","icc")
colnames(n.notes.obs) <- c("n.notes.samples","icc")
n.notes.dens.df <- rbind.data.frame(n.notes.site,n.notes.group,n.notes.obs)
head(n.notes.dens.df)
male.n.notes.plot <- ggplot(n.notes.dens.df, aes(x=n.notes.samples, fill=icc)) + geom_density()+xlim(0,1)+ #ylim(0,90)+
scale_fill_manual(name = "ICC",
values = c(alpha("blue", .5),
alpha("#6DCD59FF", .5),
alpha("#440154FF", .5)))+ theme_classic() +
guides(fill = guide_legend(size=10, keywidth = 3, title="",keyheight = 2))+ xlab("")+ylab("") +theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
axis.text.x = element_text(size=18))+ theme(legend.text=element_text(size=12,face = "bold"))
male.n.notes.plot
## Note 2 duration
note.rate.site <- cbind.data.frame(icc.site$note.rate, rep("Site",length(icc.site$note.rate)))
note.rate.group <- cbind.data.frame(icc.group$note.rate,rep("Group",length(icc.group$note.rate)))
note.rate.obs <- cbind.data.frame(icc.obs$note.rate,rep("Obs",length(icc.group$note.rate)))
colnames(note.rate.site) <- c("note.rate.samples","icc")
colnames(note.rate.group) <- c("note.rate.samples","icc")
colnames(note.rate.obs) <- c("note.rate.samples","icc")
note.rate.dens.df <- rbind.data.frame(note.rate.site,note.rate.group,note.rate.obs)
head(note.rate.dens.df)
male.note.rate.plot <- ggplot(note.rate.dens.df, aes(x=note.rate.samples, fill=icc)) + geom_density()+xlim(0,1)+ #ylim(0,90)+
scale_fill_manual(name = "ICC",
values = c(alpha("blue", .5),
alpha("#6DCD59FF", .5),
alpha("#440154FF", .5)))+ theme_classic() +
guides(fill=F)+xlab("")+ylab("") +theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
axis.text.x = element_text(size=18))
male.note.rate.plot
# Combine all plots
male.plot <- cowplot::plot_grid(male.note.1.dur.plot, male.note.2.dur.plot,male.note.1.max.freq.plot,male.note.2.max.freq.plot,
male.term.f.plot,male.n.notes.plot,male.note.rate.plot,
label_x=0.2,label_y = 1,
labels=c(" Note 1 Dur", "Note 2 Dur", "Note 1 Max Freq",
"Note 2 Max Freq","Terminal Freq","Number of Notes", " Note Rate"),
ncol=2)
cowplot::plot_grid(female.plot,male.plot,nrow=2)
# Extract the variance/co-variance matrix for each level of analysis
Vcov.site <- extract(m.stan.tarsier.f.seven.features,
pars= c("Vcov_site"),
permuted=TRUE)$Vcov_site
Vcov.female <- extract(m.stan.tarsier.f.seven.features,
pars= c("Vcov_group"),
permuted=TRUE)$Vcov_group
Vcov.obs <- extract(m.stan.tarsier.f.seven.features,
pars= c("Vcov_obs"),
permuted=TRUE)$Vcov_obs
### Turn each variance/co-variance matrix into correlation matrix
Vcov.site <- apply(X=Vcov.site, MAR=1, FUN="cov2cor")
postmean.vcov.cor.site <- apply(X=Vcov.site, MAR=1, FUN ="mean")
postmean.vcov.cor.site <- matrix(postmean.vcov.cor.site, nrow= 7, ncol=7, byrow=T)
row.names(postmean.vcov.cor.site) <-c("Note 1 Dur", "Note 2 Dur", "Note 1 Max Freq",
"Note 2 Max Freq","Terminal Freq","Number of Notes", "Note Rate")
colnames(postmean.vcov.cor.site) <- row.names(postmean.vcov.cor.site)
Vcov.female <- apply(X=Vcov.female, MAR=1, FUN="cov2cor")
postmean.vcov.cor.female <- apply(X=Vcov.female, MAR=1, FUN ="mean")
postmean.vcov.cor.female <- matrix(postmean.vcov.cor.female, nrow= 7, ncol=7, byrow=T)
row.names(postmean.vcov.cor.female) <-c("note.1.dur", "note.2.dur","note.1.max.freq",
"note.2.max.freq","term.f","n.notes","note.rate")
colnames(postmean.vcov.cor.female) <- row.names(postmean.vcov.cor.female)
Vcov.cor.obs <- apply(X=Vcov.obs, MAR=1, FUN="cov2cor")
postmean.vcov.cor.obs <- apply(X=Vcov.cor.obs, MAR=1, FUN ="mean")
postmean.vcov.cor.obs <- matrix(postmean.vcov.cor.obs, nrow= 7, ncol=7, byrow=T)
row.names(postmean.vcov.cor.obs) <- c("Note 1 Dur", "Note 2 Dur", "Note 1 Max Freq",
"Note 2 Max Freq","Terminal Freq","Number of Notes", "Note Rate")
colnames(postmean.vcov.cor.obs) <- row.names(postmean.vcov.cor.obs)
library(corrplot)
library(viridis)
## Call-level matrix
colnames(postmean.vcov.cor.obs)
corrplot(postmean.vcov.cor.obs, method=c("circle"),
col=viridis(10), cl.pos = "n",
type="lower",
#title= "Call",
tl.cex = 1.2,
tl.srt = 45, tl.col="black", diag = F, font=2)
title(main="", cex.main=2)
## Female-level matrix
corrplot(postmean.vcov.cor.female, method=c("circle"),
col=viridis(10), cl.pos = "n",
type="lower",
#title= "Site",
tl.cex = 1.2,
tl.srt = 45, tl.col="black", diag = F, font=2)
title(main="", cex.main=2)
## Site-level matrix
corrplot(postmean.vcov.cor.site,method=c("circle"),
col=viridis(10), #cl.pos = "n",
type="lower",
#title= "Site",
tl.cex = 1.2,
tl.srt = 45, tl.col="black", diag = F, font=2)
title(main="", cex.main=2)
# Extract the variance/co-variance matrix for each level of analysis
Vcov.site <- extract(m.stan.tarsier.m,
pars= c("Vcov_site"),
permuted=TRUE)$Vcov_site
Vcov.female <- extract(m.stan.tarsier.m,
pars= c("Vcov_group"),
permuted=TRUE)$Vcov_group
Vcov.obs <- extract(m.stan.tarsier.m,
pars= c("Vcov_obs"),
permuted=TRUE)$Vcov_obs
### Turn each variance/co-variance matrix into correlation matrix
Vcov.site <- apply(X=Vcov.site, MAR=1, FUN="cov2cor")
postmean.vcov.cor.site <- apply(X=Vcov.site, MAR=1, FUN ="mean")
postmean.vcov.cor.site <- matrix(postmean.vcov.cor.site, nrow= 7, ncol=7, byrow=T)
row.names(postmean.vcov.cor.site) <-c("Note 1 Dur", "Note 2 Dur", "Note 1 Max Freq",
"Note 2 Max Freq","Terminal Freq","Number of Notes", "Note Rate")
colnames(postmean.vcov.cor.site) <- row.names(postmean.vcov.cor.site)
Vcov.female <- apply(X=Vcov.female, MAR=1, FUN="cov2cor")
postmean.vcov.cor.female <- apply(X=Vcov.female, MAR=1, FUN ="mean")
postmean.vcov.cor.female <- matrix(postmean.vcov.cor.female, nrow= 7, ncol=7, byrow=T)
row.names(postmean.vcov.cor.female) <-c("note.1.dur", "note.2.dur","note.1.max.freq",
"note.2.max.freq","term.f","n.notes","note.rate")
colnames(postmean.vcov.cor.female) <- row.names(postmean.vcov.cor.female)
Vcov.cor.obs <- apply(X=Vcov.obs, MAR=1, FUN="cov2cor")
postmean.vcov.cor.obs <- apply(X=Vcov.cor.obs, MAR=1, FUN ="mean")
postmean.vcov.cor.obs <- matrix(postmean.vcov.cor.obs, nrow= 7, ncol=7, byrow=T)
row.names(postmean.vcov.cor.obs) <- c("note.1.dur", "note.2.dur","note.1.max.freq",
"note.2.max.freq","term.f","n.notes","note.rate")
colnames(postmean.vcov.cor.obs) <- row.names(postmean.vcov.cor.obs)
library(corrplot)
library(viridis)
## Call-level matrix
colnames(postmean.vcov.cor.obs)
corrplot(postmean.vcov.cor.obs, method=c("circle"),
col=viridis(10), cl.pos = "n",
type="lower",
#title= "Call",
tl.cex = 1.2,
tl.srt = 45, tl.col="black", diag = F, font=2)
title(main="", cex.main=2)
## Female-level matrix
corrplot(postmean.vcov.cor.female, method=c("circle"),
col=viridis(10), cl.pos = "n",
type="lower",
#title= "Site",
tl.cex = 1.2,
tl.srt = 45, tl.col="black", diag = F, font=2)
title(main="", cex.main=2)
## Site-level matrix
corrplot(postmean.vcov.cor.site,method=c("circle"),
col=viridis(10), #cl.pos = "n",
type="lower",
#title= "Site",
tl.cex = 1.2,
tl.srt = 45, tl.col="black", diag = F, font=2)
title(main="", cex.main=2)
DenaJGibbon/geographic_tarsiers documentation built on May 23, 2019, 12:50 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# Tarsier MANOVA
### Load required libraries
library(corrplot)
library(bayesplot)
library(rstan)
rstan_options(auto_write=TRUE)
options(mc.cores = parallel::detectCores())
library(stringr)
library(viridis)
library(ggplot2)
library(stringr)
# Read in the data
tarsier.feature.table.f <- read.csv("/Users/denasmacbook/Geographic tarsiers/tarsier.female.march28.csv")
# Isolate features of interest
d.manova <- tarsier.feature.table.f[,c("note.1.dur", "note.2.dur","note.1.max.freq",
"note.2.max.freq","term.f","n.notes","note.rate")]
## Check the structure of the data
str(d.manova)
pairs(log(d.manova))
cor(d.manova)
## Log transform data
d.manova <- log(d.manova)
hist(d.manova$note.1.dur)
### Set-up data to pass to Stan.
# Integer-coded vector of group IDs
group.int <- as.numeric(tarsier.feature.table.f$group)
# Check structure
table(group.int)
# Integer-coded vector of site IDs
site.int <- as.numeric(as.factor(tarsier.feature.table.f$site))
# Check structure
table(site.int)
# Center data matrix at feature means
col.means <- apply(d.manova, MARGIN=2, FUN="mean")
y.centered <- sweep(d.manova, MARGIN=2, STATS=col.means)
# Create a data list to pass to Stan
data_list <- list(
K = dim(d.manova)[2],
J= length(unique(tarsier.feature.table.f$group)),
M= length(unique(tarsier.feature.table.f$site)),
N= dim(d.manova)[1],
y= as.matrix(y.centered), ## features centered at zero
group= group.int,
site= site.int
)
# Code to run the STAN mdoel
# NOTE: the .stan file must be linked in the code below
m.stan.tarsier.f.seven.features = stan(file="/Users/denasmacbook/Geographic tarsiers/tarsier.stan",
model_name = "m.stan.tarsier.f",
data=data_list, iter=1500, warmup=1000, chains=2,
cores=2,
control = list(stepsize = 0.5, adapt_delta = 0.99, max_treedepth = 20))
#save(m.stan.tarsier.f, file = "/Users/denasmacbook/Desktop/m.stan.tarsier.f.rda")
#load("/Users/denasmacbook/Desktop/m.stan.tarsier.f.rda")
## Check model output
# Create traceplots to check for mixing for site level variance
draws <- as.array(m.stan.tarsier.f.seven.features, pars="ICC_site")
mcmc_trace(draws)
stan_dens(m.stan.tarsier.f.seven.features, pars=c("ICC_site"))
round(summary(m.stan.tarsier.f.seven.features, pars=c("ICC_site"))$summary, 3)
# Create traceplots to check for mixing for group level variance
draws <- as.array(m.stan.tarsier.f.seven.features, pars="ICC_group")
mcmc_trace(draws)
stan_dens(m.stan.tarsier.f.seven.features, pars=c("ICC_group"))
round(summary(m.stan.tarsier.f.seven.features, pars=c("ICC_group"))$summary, 3)
# Check degrees of freedom parameter
stan_dens(m.stan.tarsier.f.seven.features, pars=c("DF_obs"))
round(summary(m.stan.tarsier.f.seven.features, pars=c("DF_obs"))$summary, 3)
stan_dens(m.stan.tarsier.f, pars=c("DF_group"))
round(summary(m.stan.tarsier.f, pars=c("DF_group"))$summary, 3)
stan_dens(m.stan.tarsier.f, pars=c("DF_site"))
round(summary(m.stan.tarsier.f, pars=c("DF_site"))$summary, 3)
post.samples <- rstan::extract(m.stan.tarsier.f.seven.features,
pars= c("ICC_site", "ICC_group", "Maha_sqd", "DF_site","DF_group",
"DF_obs", "Vcov_site", "Vcov_group", "Vcov_obs"),
permuted=TRUE)
# Create vector with feature names
features <- c("note.1.dur", "note.2.dur","note.1.max.freq",
"note.2.max.freq","term.f","n.notes","note.rate")
# Add column names to samples
colnames(post.samples$ICC_group) <- features
colnames(post.samples$ICC_site) <- features
# Convert to dataframe
icc.site <- as.data.frame(post.samples$ICC_site)
icc.group <- as.data.frame(post.samples$ICC_group)
# Calculate ICC for observation-level and add column names
icc.obs <- as.data.frame(1 - icc.group - icc.site)
colnames(icc.obs) <- features
### Code to make Figure 4. Posterior densities for the intraclass correlation coefficients for the three levels in our dataset
### (call, fefemale and site) for each feature of the tarsier duet phrase
## Note 2 duration
note.1.dur.site <- cbind.data.frame(icc.site$note.1.dur, rep("Site",length(icc.site$note.1.dur)))
note.1.dur.group <- cbind.data.frame(icc.group$note.1.dur,rep("Group",length(icc.group$note.1.dur)))
note.1.dur.obs <- cbind.data.frame(icc.obs$note.1.dur,rep("Obs",length(icc.obs$note.1.dur)))
colnames(note.1.dur.site) <- c("note.1.dur","icc")
colnames(note.1.dur.group) <- c("note.1.dur","icc")
colnames(note.1.dur.obs) <- c("note.1.dur","icc")
note.1.dur.dens.df <- rbind.data.frame(note.1.dur.site,note.1.dur.group,note.1.dur.obs)
head(note.1.dur.dens.df)
female.note.1.dur.plot <- ggplot(note.1.dur.dens.df, aes(x=note.1.dur, fill=icc)) + geom_density()+xlim(0,1)+#+ylim(0,4)+
scale_fill_manual(name = "icc",
values = c(alpha("blue", .5),
alpha("#6DCD59FF", .5),
alpha("#440154FF", .5)))+ theme_classic() +
guides(fill=F)+xlab("")+ylab("") +theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
axis.text.x = element_text(size=18))
female.note.1.dur.plot
## Note 2 duration
note.2.dur.site <- cbind.data.frame(icc.site$note.2.dur, rep("Site",length(icc.site$note.2.dur)))
note.2.dur.group <- cbind.data.frame(icc.group$note.2.dur,rep("Group",length(icc.group$note.2.dur)))
note.2.dur.obs <- cbind.data.frame(icc.obs$note.2.dur,rep("Obs",length(icc.obs$note.2.dur)))
colnames(note.2.dur.site) <- c("note.2.dur","icc")
colnames(note.2.dur.group) <- c("note.2.dur","icc")
colnames(note.2.dur.obs) <- c("note.2.dur","icc")
note.2.dur.dens.df <- rbind.data.frame(note.2.dur.site,note.2.dur.group,note.2.dur.obs)
head(note.2.dur.dens.df)
female.note.2.dur.plot <- ggplot(note.2.dur.dens.df, aes(x=note.2.dur, fill=icc)) + geom_density()+xlim(0,1)+#+ylim(0,4)+
scale_fill_manual(name = "icc",
values = c(alpha("blue", .5),
alpha("#6DCD59FF", .5),
alpha("#440154FF", .5)))+ theme_classic() +
guides(fill=F)+xlab("")+ylab("") +theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
axis.text.x = element_text(size=18))
female.note.2.dur.plot
## Note 1 Max Freq
note.1.max.freq.site <- cbind.data.frame(icc.site$note.1.max.freq, rep("Site",length(icc.site$note.1.max.freq)))
note.1.max.freq.group <- cbind.data.frame(icc.group$note.1.max.freq,rep("Group",length(icc.group$note.1.max.freq)))
note.1.max.freq.obs <- cbind.data.frame(icc.obs$note.1.max.freq,rep("Obs",length(icc.obs$note.1.max.freq)))
colnames(note.1.max.freq.site) <- c("note.1.max.freq.samples","icc")
colnames(note.1.max.freq.group) <- c("note.1.max.freq.samples","icc")
colnames(note.1.max.freq.obs) <- c("note.1.max.freq.samples","icc")
note.1.max.freq.dens.df <- rbind.data.frame(note.1.max.freq.site,note.1.max.freq.group,note.1.max.freq.obs)
head(note.1.max.freq.dens.df)
female.note.1.max.freq.plot <- ggplot(note.1.max.freq.dens.df, aes(x=note.1.max.freq.samples, fill=icc)) + geom_density()+xlim(0,1)+#+ylim(0,4)+
scale_fill_manual(name = "icc",
values = c(alpha("blue", .5),
alpha("#6DCD59FF", .5),
alpha("#440154FF", .5)))+ theme_classic() +
guides(fill=F)+xlab("")+ylab("") +theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
axis.text.x = element_text(size=18))
female.note.1.max.freq.plot
# Note 2 Max Freq
note.2.max.freq.site <- cbind.data.frame(icc.site$note.2.max.freq, rep("Site",length(icc.site$note.2.max.freq)))
note.2.max.freq.group <- cbind.data.frame(icc.group$note.2.max.freq,rep("Group",length(icc.group$note.2.max.freq)))
note.2.max.freq.obs <- cbind.data.frame(icc.obs$note.2.max.freq,rep("Obs",length(icc.obs$note.2.max.freq)))
colnames(note.2.max.freq.site) <- c("note.2.max.freq.samples","icc")
colnames(note.2.max.freq.group) <- c("note.2.max.freq.samples","icc")
colnames(note.2.max.freq.obs) <- c("note.2.max.freq.samples","icc")
note.2.max.freq.dens.df <- rbind.data.frame(note.2.max.freq.site,note.2.max.freq.group,note.2.max.freq.obs)
head(note.2.max.freq.dens.df)
female.note.2.max.freq.plot <- ggplot(note.2.max.freq.dens.df, aes(x=note.2.max.freq.samples, fill=icc)) + geom_density()+xlim(0,1)+ylim(0,4)+
scale_fill_manual(name = "icc",
values = c(alpha("blue", .5),
alpha("#6DCD59FF", .5),
alpha("#440154FF", .5)))+ theme_classic() +
guides(fill=F)+xlab("")+ylab("") +theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
axis.text.x = element_text(size=18))
female.note.2.max.freq.plot
## Note 1 max freq
term.f.site <- cbind.data.frame(icc.site$term.f, rep("Between-site",length(icc.site$term.f)))
term.f.group <- cbind.data.frame(icc.group$term.f,rep("Between-fefemale",length(icc.group$term.f)))
term.f.obs <- cbind.data.frame(icc.obs$term.f,rep("Within-fefemale",length(icc.group$term.f)))
colnames(term.f.site) <- c("term.f.samples","icc")
colnames(term.f.group) <- c("term.f.samples","icc")
colnames(term.f.obs) <- c("term.f.samples","icc")
term.f.dens.df <- rbind.data.frame(term.f.site,term.f.group,term.f.obs)
head(term.f.dens.df)
female.term.f.plot <- ggplot(term.f.dens.df, aes(x=term.f.samples, fill=icc)) + geom_density()+xlim(0,1)+ #ylim(0,90)+
scale_fill_manual(name = "ICC",
values = c(alpha("blue", .5),
alpha("#6DCD59FF", .5),
alpha("#440154FF", .5)))+ theme_classic() +
guides(fill=F)+xlab("")+ylab("") +theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
axis.text.x = element_text(size=18),legend.text=element_text(size=24))+
xlab("")+ylab("") + theme(legend.position = c(.8,.87))
female.term.f.plot
## Note 2 duration
n.notes.site <- cbind.data.frame(icc.site$n.notes, rep("Island",length(icc.site$n.notes)))
n.notes.group <- cbind.data.frame(icc.group$n.notes,rep("Female",length(icc.group$n.notes)))
n.notes.obs <- cbind.data.frame(icc.obs$n.notes,rep("Phrase",length(icc.group$n.notes)))
colnames(n.notes.site) <- c("n.notes.samples","icc")
colnames(n.notes.group) <- c("n.notes.samples","icc")
colnames(n.notes.obs) <- c("n.notes.samples","icc")
n.notes.dens.df <- rbind.data.frame(n.notes.site,n.notes.group,n.notes.obs)
head(n.notes.dens.df)
female.n.notes.plot <- ggplot(n.notes.dens.df, aes(x=n.notes.samples, fill=icc)) + geom_density()+xlim(0,1)+ #ylim(0,90)+
scale_fill_manual(name = "ICC",
values = c(alpha("blue", .5),
alpha("#6DCD59FF", .5),
alpha("#440154FF", .5)))+ theme_classic() +
guides(fill = guide_legend(size=10, keywidth = 3, title="",keyheight = 2))+ xlab("")+ylab("") +theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
axis.text.x = element_text(size=18))+ theme(legend.text=element_text(size=12,face = "bold"))
female.n.notes.plot
## Note 2 duration
note.rate.site <- cbind.data.frame(icc.site$note.rate, rep("Site",length(icc.site$note.rate)))
note.rate.group <- cbind.data.frame(icc.group$note.rate,rep("Group",length(icc.group$note.rate)))
note.rate.obs <- cbind.data.frame(icc.obs$note.rate,rep("Obs",length(icc.group$note.rate)))
colnames(note.rate.site) <- c("note.rate.samples","icc")
colnames(note.rate.group) <- c("note.rate.samples","icc")
colnames(note.rate.obs) <- c("note.rate.samples","icc")
note.rate.dens.df <- rbind.data.frame(note.rate.site,note.rate.group,note.rate.obs)
head(note.rate.dens.df)
female.note.rate.plot <- ggplot(note.rate.dens.df, aes(x=note.rate.samples, fill=icc)) + geom_density()+xlim(0,1)+ #ylim(0,90)+
scale_fill_manual(name = "ICC",
values = c(alpha("blue", .5),
alpha("#6DCD59FF", .5),
alpha("#440154FF", .5)))+ theme_classic() +
guides(fill=F)+xlab("")+ylab("") +theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
axis.text.x = element_text(size=18))
female.note.rate.plot
# Combine all plots
female.plot <- cowplot::plot_grid(female.note.1.dur.plot, female.note.2.dur.plot,female.note.1.max.freq.plot,female.note.2.max.freq.plot,
female.term.f.plot,female.n.notes.plot,female.note.rate.plot,
label_x=0.3,label_y = 1,
labels=c(" Note 1 Dur", "Note 2 Dur", "Note 1 Max Freq",
"Note 2 Max Freq","Terminal Freq","Number of Notes", " Note Rate"),
ncol=2)
######## Male MANOVA
male.tarsier.df <- read.csv("tarsier.feature.table.m.csv")
d.manova.male <- male.tarsier.df[,c("note.1.dur", "note.2.dur","note.1.max.freq",
"note.2.max.freq","term.f","n.notes","note.rate")]
## Check the structure of the data
str(d.manova.male)
pairs(d.manova.male)
cor(d.manova.male)
## Log transform data
d.manova.male <- log(d.manova.male)
### Set-up data to pass to Stan.
# Integer-coded vector of group IDs
group.int <- as.numeric(male.tarsier.df$group)
# Check structure
table(group.int)
# Integer-coded vector of site IDs
site.int <- as.numeric(as.factor(male.tarsier.df$site))
# Check structure
table(site.int)
# Center data matrix at feature means
col.means <- apply(d.manova.male, MARGIN=2, FUN="mean")
y.centered <- sweep(d.manova.male, MARGIN=2, STATS=col.means)
# Create a data list to pass to Stan
data_list <- list(
K = dim(d.manova.male)[2],
J= length(unique(male.tarsier.df$group)),
M= length(unique(male.tarsier.df$site)),
N= dim(d.manova.male)[1],
y= as.matrix(y.centered), ## features centered at zero
group= group.int,
site= site.int
)
# Code to run the STAN mdoel
# NOTE: the .stan file must be linked in the code below
m.stan.tarsier.m = stan(file="/Users/denasmacbook/Desktop/Clink et al Multivariate Analysis/MANOVA_MNG_oct18_2017.stan",
model_name = "m.stan.tarsier.m",
data=data_list, iter=1000, warmup=500, chains=2,
cores=2,
control = list(stepsize = 0.5, adapt_delta = 0.99, max_treedepth = 20))
## Check model output
# Create traceplots to check for mixing for site level variance
draws <- as.array(m.stan.tarsier.m, pars="ICC_site")
mcmc_trace(draws)
stan_dens(m.stan.tarsier.m, pars=c("ICC_site"))
round(summary(m.stan.tarsier.m, pars=c("ICC_site"))$summary, 3)
# Create traceplots to check for mixing for group level variance
draws <- as.array(m.stan.tarsier.m, pars="ICC_group")
mcmc_trace(draws)
stan_dens(m.stan.tarsier.m, pars=c("ICC_group"))
round(summary(m.stan.tarsier.m, pars=c("ICC_group"))$summary, 3)
# Check degrees of freedom parameter
stan_dens(m.stan.tarsier.m, pars=c("DF_obs"))
round(summary(m.stan.tarsier.m, pars=c("DF_obs"))$summary, 3)
stan_dens(m.stan.tarsier.m, pars=c("DF_group"))
round(summary(m.stan.tarsier.m, pars=c("DF_group"))$summary, 3)
stan_dens(m.stan.tarsier.m, pars=c("DF_site"))
round(summary(m.stan.tarsier.m, pars=c("DF_site"))$summary, 3)
post.samples <- extract(m.stan.tarsier.m,
pars= c("ICC_site", "ICC_group", "Maha_sqd", "DF_site","DF_group",
"DF_obs", "Vcov_site", "Vcov_group", "Vcov_obs"),
permuted=TRUE)
# Create vector with feature names
features <- c("note.1.dur", "note.2.dur","note.1.max.freq",
"note.2.max.freq","term.f","n.notes","note.rate")
# Add column names to samples
colnames(post.samples$ICC_group) <- features
colnames(post.samples$ICC_site) <- features
# Convert to dataframe
icc.site <- as.data.frame(post.samples$ICC_site)
icc.group <- as.data.frame(post.samples$ICC_group)
# Calculate ICC for observation-level and add column names
icc.obs <- as.data.frame(1 - icc.group - icc.site)
colnames(icc.obs) <- features
### Code to make Figure 4. Posterior densities for the intraclass correlation coefficients for the three levels in our dataset
### (call, female and site) for each feature of the tarsier duet phrase
## Note 2 duration
note.1.dur.site <- cbind.data.frame(icc.site$note.1.dur, rep("Site",length(icc.site$note.1.dur)))
note.1.dur.group <- cbind.data.frame(icc.group$note.1.dur,rep("Group",length(icc.group$note.1.dur)))
note.1.dur.obs <- cbind.data.frame(icc.obs$note.1.dur,rep("Obs",length(icc.obs$note.1.dur)))
colnames(note.1.dur.site) <- c("note.1.dur","icc")
colnames(note.1.dur.group) <- c("note.1.dur","icc")
colnames(note.1.dur.obs) <- c("note.1.dur","icc")
note.1.dur.dens.df <- rbind.data.frame(note.1.dur.site,note.1.dur.group,note.1.dur.obs)
head(note.1.dur.dens.df)
male.note.1.dur.plot <- ggplot(note.1.dur.dens.df, aes(x=note.1.dur, fill=icc)) + geom_density()+xlim(0,1)+#+ylim(0,4)+
scale_fill_manual(name = "icc",
values = c(alpha("blue", .5),
alpha("#6DCD59FF", .5),
alpha("#440154FF", .5)))+ theme_classic() +
guides(fill=F)+xlab("")+ylab("") +theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
axis.text.x = element_text(size=18))
male.note.1.dur.plot
## Note 2 duration
note.2.dur.site <- cbind.data.frame(icc.site$note.2.dur, rep("Site",length(icc.site$note.2.dur)))
note.2.dur.group <- cbind.data.frame(icc.group$note.2.dur,rep("Group",length(icc.group$note.2.dur)))
note.2.dur.obs <- cbind.data.frame(icc.obs$note.2.dur,rep("Obs",length(icc.obs$note.2.dur)))
colnames(note.2.dur.site) <- c("note.2.dur","icc")
colnames(note.2.dur.group) <- c("note.2.dur","icc")
colnames(note.2.dur.obs) <- c("note.2.dur","icc")
note.2.dur.dens.df <- rbind.data.frame(note.2.dur.site,note.2.dur.group,note.2.dur.obs)
head(note.2.dur.dens.df)
male.note.2.dur.plot <- ggplot(note.2.dur.dens.df, aes(x=note.2.dur, fill=icc)) + geom_density()+xlim(0,1)+#+ylim(0,4)+
scale_fill_manual(name = "icc",
values = c(alpha("blue", .5),
alpha("#6DCD59FF", .5),
alpha("#440154FF", .5)))+ theme_classic() +
guides(fill=F)+xlab("")+ylab("") +theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
axis.text.x = element_text(size=18))
male.note.2.dur.plot
## Note 1 Max Freq
note.1.max.freq.site <- cbind.data.frame(icc.site$note.1.max.freq, rep("Site",length(icc.site$note.1.max.freq)))
note.1.max.freq.group <- cbind.data.frame(icc.group$note.1.max.freq,rep("Group",length(icc.group$note.1.max.freq)))
note.1.max.freq.obs <- cbind.data.frame(icc.obs$note.1.max.freq,rep("Obs",length(icc.obs$note.1.max.freq)))
colnames(note.1.max.freq.site) <- c("note.1.max.freq.samples","icc")
colnames(note.1.max.freq.group) <- c("note.1.max.freq.samples","icc")
colnames(note.1.max.freq.obs) <- c("note.1.max.freq.samples","icc")
note.1.max.freq.dens.df <- rbind.data.frame(note.1.max.freq.site,note.1.max.freq.group,note.1.max.freq.obs)
head(note.1.max.freq.dens.df)
male.note.1.max.freq.plot <- ggplot(note.1.max.freq.dens.df, aes(x=note.1.max.freq.samples, fill=icc)) + geom_density()+xlim(0,1)+#+ylim(0,4)+
scale_fill_manual(name = "icc",
values = c(alpha("blue", .5),
alpha("#6DCD59FF", .5),
alpha("#440154FF", .5)))+ theme_classic() +
guides(fill=F)+xlab("")+ylab("") +theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
axis.text.x = element_text(size=18))
male.note.1.max.freq.plot
# Note 2 Max Freq
note.2.max.freq.site <- cbind.data.frame(icc.site$note.2.max.freq, rep("Site",length(icc.site$note.2.max.freq)))
note.2.max.freq.group <- cbind.data.frame(icc.group$note.2.max.freq,rep("Group",length(icc.group$note.2.max.freq)))
note.2.max.freq.obs <- cbind.data.frame(icc.obs$note.2.max.freq,rep("Obs",length(icc.obs$note.2.max.freq)))
colnames(note.2.max.freq.site) <- c("note.2.max.freq.samples","icc")
colnames(note.2.max.freq.group) <- c("note.2.max.freq.samples","icc")
colnames(note.2.max.freq.obs) <- c("note.2.max.freq.samples","icc")
note.2.max.freq.dens.df <- rbind.data.frame(note.2.max.freq.site,note.2.max.freq.group,note.2.max.freq.obs)
head(note.2.max.freq.dens.df)
male.note.2.max.freq.plot <- ggplot(note.2.max.freq.dens.df, aes(x=note.2.max.freq.samples, fill=icc)) + geom_density()+xlim(0,1)+#+ylim(0,4)+
scale_fill_manual(name = "icc",
values = c(alpha("blue", .5),
alpha("#6DCD59FF", .5),
alpha("#440154FF", .5)))+ theme_classic() +
guides(fill=F)+xlab("")+ylab("") +theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
axis.text.x = element_text(size=18))
male.note.2.max.freq.plot
## Note 1 max freq
term.f.site <- cbind.data.frame(icc.site$term.f, rep("Between-site",length(icc.site$term.f)))
term.f.group <- cbind.data.frame(icc.group$term.f,rep("Between-female",length(icc.group$term.f)))
term.f.obs <- cbind.data.frame(icc.obs$term.f,rep("Within-female",length(icc.group$term.f)))
colnames(term.f.site) <- c("term.f.samples","icc")
colnames(term.f.group) <- c("term.f.samples","icc")
colnames(term.f.obs) <- c("term.f.samples","icc")
term.f.dens.df <- rbind.data.frame(term.f.site,term.f.group,term.f.obs)
head(term.f.dens.df)
male.term.f.plot <- ggplot(term.f.dens.df, aes(x=term.f.samples, fill=icc)) + geom_density()+xlim(0,1)+ #ylim(0,90)+
scale_fill_manual(name = "ICC",
values = c(alpha("blue", .5),
alpha("#6DCD59FF", .5),
alpha("#440154FF", .5)))+ theme_classic() +
guides(fill=F)+xlab("")+ylab("") +theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
axis.text.x = element_text(size=18),legend.text=element_text(size=24))+
xlab("")+ylab("") + theme(legend.position = c(.8,.87))
male.term.f.plot
## Note 2 duration
n.notes.site <- cbind.data.frame(icc.site$n.notes, rep("Island",length(icc.site$n.notes)))
n.notes.group <- cbind.data.frame(icc.group$n.notes,rep("Male",length(icc.group$n.notes)))
n.notes.obs <- cbind.data.frame(icc.obs$n.notes,rep("Phrase",length(icc.group$n.notes)))
colnames(n.notes.site) <- c("n.notes.samples","icc")
colnames(n.notes.group) <- c("n.notes.samples","icc")
colnames(n.notes.obs) <- c("n.notes.samples","icc")
n.notes.dens.df <- rbind.data.frame(n.notes.site,n.notes.group,n.notes.obs)
head(n.notes.dens.df)
male.n.notes.plot <- ggplot(n.notes.dens.df, aes(x=n.notes.samples, fill=icc)) + geom_density()+xlim(0,1)+ #ylim(0,90)+
scale_fill_manual(name = "ICC",
values = c(alpha("blue", .5),
alpha("#6DCD59FF", .5),
alpha("#440154FF", .5)))+ theme_classic() +
guides(fill = guide_legend(size=10, keywidth = 3, title="",keyheight = 2))+ xlab("")+ylab("") +theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
axis.text.x = element_text(size=18))+ theme(legend.text=element_text(size=12,face = "bold"))
male.n.notes.plot
## Note 2 duration
note.rate.site <- cbind.data.frame(icc.site$note.rate, rep("Site",length(icc.site$note.rate)))
note.rate.group <- cbind.data.frame(icc.group$note.rate,rep("Group",length(icc.group$note.rate)))
note.rate.obs <- cbind.data.frame(icc.obs$note.rate,rep("Obs",length(icc.group$note.rate)))
colnames(note.rate.site) <- c("note.rate.samples","icc")
colnames(note.rate.group) <- c("note.rate.samples","icc")
colnames(note.rate.obs) <- c("note.rate.samples","icc")
note.rate.dens.df <- rbind.data.frame(note.rate.site,note.rate.group,note.rate.obs)
head(note.rate.dens.df)
male.note.rate.plot <- ggplot(note.rate.dens.df, aes(x=note.rate.samples, fill=icc)) + geom_density()+xlim(0,1)+ #ylim(0,90)+
scale_fill_manual(name = "ICC",
values = c(alpha("blue", .5),
alpha("#6DCD59FF", .5),
alpha("#440154FF", .5)))+ theme_classic() +
guides(fill=F)+xlab("")+ylab("") +theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
axis.text.x = element_text(size=18))
male.note.rate.plot
# Combine all plots
male.plot <- cowplot::plot_grid(male.note.1.dur.plot, male.note.2.dur.plot,male.note.1.max.freq.plot,male.note.2.max.freq.plot,
male.term.f.plot,male.n.notes.plot,male.note.rate.plot,
label_x=0.2,label_y = 1,
labels=c(" Note 1 Dur", "Note 2 Dur", "Note 1 Max Freq",
"Note 2 Max Freq","Terminal Freq","Number of Notes", " Note Rate"),
ncol=2)
cowplot::plot_grid(female.plot,male.plot,nrow=2)
# Extract the variance/co-variance matrix for each level of analysis
Vcov.site <- extract(m.stan.tarsier.f.seven.features,
pars= c("Vcov_site"),
permuted=TRUE)$Vcov_site
Vcov.female <- extract(m.stan.tarsier.f.seven.features,
pars= c("Vcov_group"),
permuted=TRUE)$Vcov_group
Vcov.obs <- extract(m.stan.tarsier.f.seven.features,
pars= c("Vcov_obs"),
permuted=TRUE)$Vcov_obs
### Turn each variance/co-variance matrix into correlation matrix
Vcov.site <- apply(X=Vcov.site, MAR=1, FUN="cov2cor")
postmean.vcov.cor.site <- apply(X=Vcov.site, MAR=1, FUN ="mean")
postmean.vcov.cor.site <- matrix(postmean.vcov.cor.site, nrow= 7, ncol=7, byrow=T)
row.names(postmean.vcov.cor.site) <-c("Note 1 Dur", "Note 2 Dur", "Note 1 Max Freq",
"Note 2 Max Freq","Terminal Freq","Number of Notes", "Note Rate")
colnames(postmean.vcov.cor.site) <- row.names(postmean.vcov.cor.site)
Vcov.female <- apply(X=Vcov.female, MAR=1, FUN="cov2cor")
postmean.vcov.cor.female <- apply(X=Vcov.female, MAR=1, FUN ="mean")
postmean.vcov.cor.female <- matrix(postmean.vcov.cor.female, nrow= 7, ncol=7, byrow=T)
row.names(postmean.vcov.cor.female) <-c("note.1.dur", "note.2.dur","note.1.max.freq",
"note.2.max.freq","term.f","n.notes","note.rate")
colnames(postmean.vcov.cor.female) <- row.names(postmean.vcov.cor.female)
Vcov.cor.obs <- apply(X=Vcov.obs, MAR=1, FUN="cov2cor")
postmean.vcov.cor.obs <- apply(X=Vcov.cor.obs, MAR=1, FUN ="mean")
postmean.vcov.cor.obs <- matrix(postmean.vcov.cor.obs, nrow= 7, ncol=7, byrow=T)
row.names(postmean.vcov.cor.obs) <- c("Note 1 Dur", "Note 2 Dur", "Note 1 Max Freq",
"Note 2 Max Freq","Terminal Freq","Number of Notes", "Note Rate")
colnames(postmean.vcov.cor.obs) <- row.names(postmean.vcov.cor.obs)
library(corrplot)
library(viridis)
## Call-level matrix
colnames(postmean.vcov.cor.obs)
corrplot(postmean.vcov.cor.obs, method=c("circle"),
col=viridis(10), cl.pos = "n",
type="lower",
#title= "Call",
tl.cex = 1.2,
tl.srt = 45, tl.col="black", diag = F, font=2)
title(main="", cex.main=2)
## Female-level matrix
corrplot(postmean.vcov.cor.female, method=c("circle"),
col=viridis(10), cl.pos = "n",
type="lower",
#title= "Site",
tl.cex = 1.2,
tl.srt = 45, tl.col="black", diag = F, font=2)
title(main="", cex.main=2)
## Site-level matrix
corrplot(postmean.vcov.cor.site,method=c("circle"),
col=viridis(10), #cl.pos = "n",
type="lower",
#title= "Site",
tl.cex = 1.2,
tl.srt = 45, tl.col="black", diag = F, font=2)
title(main="", cex.main=2)
# Extract the variance/co-variance matrix for each level of analysis
Vcov.site <- extract(m.stan.tarsier.m,
pars= c("Vcov_site"),
permuted=TRUE)$Vcov_site
Vcov.female <- extract(m.stan.tarsier.m,
pars= c("Vcov_group"),
permuted=TRUE)$Vcov_group
Vcov.obs <- extract(m.stan.tarsier.m,
pars= c("Vcov_obs"),
permuted=TRUE)$Vcov_obs
### Turn each variance/co-variance matrix into correlation matrix
Vcov.site <- apply(X=Vcov.site, MAR=1, FUN="cov2cor")
postmean.vcov.cor.site <- apply(X=Vcov.site, MAR=1, FUN ="mean")
postmean.vcov.cor.site <- matrix(postmean.vcov.cor.site, nrow= 7, ncol=7, byrow=T)
row.names(postmean.vcov.cor.site) <-c("Note 1 Dur", "Note 2 Dur", "Note 1 Max Freq",
"Note 2 Max Freq","Terminal Freq","Number of Notes", "Note Rate")
colnames(postmean.vcov.cor.site) <- row.names(postmean.vcov.cor.site)
Vcov.female <- apply(X=Vcov.female, MAR=1, FUN="cov2cor")
postmean.vcov.cor.female <- apply(X=Vcov.female, MAR=1, FUN ="mean")
postmean.vcov.cor.female <- matrix(postmean.vcov.cor.female, nrow= 7, ncol=7, byrow=T)
row.names(postmean.vcov.cor.female) <-c("note.1.dur", "note.2.dur","note.1.max.freq",
"note.2.max.freq","term.f","n.notes","note.rate")
colnames(postmean.vcov.cor.female) <- row.names(postmean.vcov.cor.female)
Vcov.cor.obs <- apply(X=Vcov.obs, MAR=1, FUN="cov2cor")
postmean.vcov.cor.obs <- apply(X=Vcov.cor.obs, MAR=1, FUN ="mean")
postmean.vcov.cor.obs <- matrix(postmean.vcov.cor.obs, nrow= 7, ncol=7, byrow=T)
row.names(postmean.vcov.cor.obs) <- c("note.1.dur", "note.2.dur","note.1.max.freq",
"note.2.max.freq","term.f","n.notes","note.rate")
colnames(postmean.vcov.cor.obs) <- row.names(postmean.vcov.cor.obs)
library(corrplot)
library(viridis)
## Call-level matrix
colnames(postmean.vcov.cor.obs)
corrplot(postmean.vcov.cor.obs, method=c("circle"),
col=viridis(10), cl.pos = "n",
type="lower",
#title= "Call",
tl.cex = 1.2,
tl.srt = 45, tl.col="black", diag = F, font=2)
title(main="", cex.main=2)
## Female-level matrix
corrplot(postmean.vcov.cor.female, method=c("circle"),
col=viridis(10), cl.pos = "n",
type="lower",
#title= "Site",
tl.cex = 1.2,
tl.srt = 45, tl.col="black", diag = F, font=2)
title(main="", cex.main=2)
## Site-level matrix
corrplot(postmean.vcov.cor.site,method=c("circle"),
col=viridis(10), #cl.pos = "n",
type="lower",
#title= "Site",
tl.cex = 1.2,
tl.srt = 45, tl.col="black", diag = F, font=2)
title(main="", cex.main=2)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.